library(caret) # load this first as it loads plyr
library(gbm)
library(car)
library(grid)
library(gridExtra)
library(randomForest)
library(xtable)
library(kableExtra)
library(abt)
library(tidyverse)
source('FISH_functions.R')
Raw fish data
fish = read.csv('data/Fish benthic physical sitelevel 2019.csv', strip.white=TRUE)
fish %>% glimpse
## Observations: 619
## Variables: 300
## $ YEAR <int> 2007, 2007, 2007, 2007, 2007, 2007, …
## $ REGION <fct> Keppel, Keppel, Keppel, Keppel, Kepp…
## $ EXPOSURE <fct> Semi-Exposed, Exposed, Semi-Exposed,…
## $ NTR <fct> Fished, Fished, Fished, Fished, Fish…
## $ NTR.Pooled <fct> Fished, Fished, Fished, Fished, Fish…
## $ SITE <fct> GK1, GK2, GK3, GK4, GK7, GK8, H3, MI…
## $ lut.arge <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ lut.carp <dbl> 13.3320, 11.3322, 4.6662, 8.6658, 16…
## $ lut.flma <dbl> 0.0000, 3.3330, 0.0000, 0.0000, 7.33…
## $ lut.fulv <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ lut.kas <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ lut.lemn <dbl> 0.0000, 0.6666, 0.0000, 0.0000, 0.00…
## $ lut.lutj <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ lut.mono <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ lut.quin <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ lut.rivu <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ lut.russ <dbl> 67.3266, 4.6662, 0.0000, 5.3328, 116…
## $ lut.seba <dbl> 1.9998, 0.6666, 9.9990, 0.0000, 0.00…
## $ lut.vitt <dbl> 377.9622, 10.6656, 16.6650, 9.3324, …
## $ mac.macu <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pagrus <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sym.nema <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ par.barb <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ par.boides <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ par.bifa <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ par.indi <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ par.multi <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ par.cili <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ upe.trag <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ platax <dbl> 0.0000, 0.6666, 0.0000, 0.6666, 0.66…
## $ Platycephalus.spp. <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ let.atki <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ let.hara <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ let.lati <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ let.lent <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ let.mini <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ let.nebu <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ let.oliv <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ let.obso <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ let.orna <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ gym.spp <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ mon.gran <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sco.bili <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sco.marg <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ sco.mono <dbl> 4.6662, 1.9998, 3.9996, 5.3328, 3.99…
## $ sco.line <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ any.leuc <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ ath.roga <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cep.argu <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cep.boen <dbl> 21.9978, 5.9994, 7.3326, 8.6658, 3.3…
## $ cep.cyan <dbl> 1.3332, 0.0000, 1.9998, 1.3332, 0.00…
## $ cep.micr <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cro.alti <dbl> 0.6666, 0.0000, 0.6666, 0.6666, 0.00…
## $ epi.caer <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ epi.coio <dbl> 0.0000, 0.0000, 0.0000, 1.3332, 0.00…
## $ epi.cora <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ epi.fasc <dbl> 0.0000, 0.0000, 0.6666, 0.0000, 0.00…
## $ epi.fusc <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ epi.howl <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ epi.lanc <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ epi.hexa <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ epi.merr <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ epi.ongu <dbl> 0.6666, 0.0000, 0.6666, 0.0000, 0.00…
## $ epi.quoy <dbl> 6.6660, 2.6664, 12.6654, 1.9998, 1.9…
## $ epi.sexf <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pms.laev <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pms.leop <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.66…
## $ pms.macu <dbl> 5.9994, 3.9996, 17.3316, 7.3326, 9.9…
## $ pte.ante <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pte.voli <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ dip.bifa <dbl> 1.9998, 1.9998, 7.9992, 1.9998, 1.99…
## $ dia.pict <dbl> 23.3310, 3.3330, 4.6662, 1.9998, 9.9…
## $ ple.chae <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ ple.flav <dbl> 1.9998, 1.3332, 2.6664, 0.6666, 0.00…
## $ ple.gibb <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ ple.less <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ ple.line <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ ple.picu <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ ple.unic <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ bodianus <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ oxy.diag <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ che.chlo <dbl> 1.3332, 0.6666, 5.9994, 0.6666, 1.33…
## $ che.fasc <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ che.tril <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ che.undu <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cho.anch <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cho.ceph <dbl> 0.6666, 0.0000, 0.6666, 1.3332, 0.66…
## $ cho.cya <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ cho.fasc <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ cho.grap <dbl> 0.0000, 0.0000, 0.6666, 0.0000, 1.33…
## $ cho.scho <dbl> 1.3332, 0.0000, 0.6666, 0.6666, 0.00…
## $ cho.vitt <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ epb.insi <dbl> 0.6666, 0.0000, 0.0000, 0.0000, 1.33…
## $ gom.vari <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ hem.melt <dbl> 17.3316, 9.3324, 25.9974, 13.9986, 2…
## $ hem.fasc <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ psa.wai <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ aca.bloc <dbl> 0.0000, 0.0000, 0.6666, 0.0000, 4.66…
## $ aca.duss <dbl> 2.6664, 0.0000, 0.0000, 0.0000, 9.99…
## $ aca.gram <dbl> 0.0000, 0.0000, 0.6666, 1.3332, 4.66…
## $ aca.line <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ aca.nans <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ aca.ncus <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ aca.nuda <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ aca.xant <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ aca.thom <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ aca.spp <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cte.bino <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cte.stri <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ nas.annu <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ nas.brach <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ nas.brevi <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ nas.litu <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ nas.tong <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ nas.unic <dbl> 0.0000, 1.9998, 0.0000, 0.0000, 0.00…
## $ pri.spp <dbl> 0.0000, 0.0000, 0.0000, 3.9996, 18.6…
## $ kyp.spp <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ zeb.scop <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ zeb.veli <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ zan.corn <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ bol.muri <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cal.caro <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cet.bico <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chl.blee <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chl.micr <dbl> 0.6666, 0.0000, 0.0000, 0.0000, 0.00…
## $ chs.sord <dbl> 10.6656, 2.6664, 0.0000, 0.0000, 0.0…
## $ hip.long <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.alti <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.cham <dbl> 6.6660, 0.6666, 0.0000, 0.0000, 0.00…
## $ sca.dimi <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.flav <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.fors <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.fren <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ sca.ghob <dbl> 7.9992, 3.9996, 1.9998, 0.6666, 13.3…
## $ sca.glob <dbl> 0.6666, 7.9992, 0.0000, 0.0000, 0.66…
## $ sca.nigr <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.ovic <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.psit <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.rivu <dbl> 37.3296, 25.3308, 19.3314, 7.9992, 4…
## $ sca.rubr <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.schl <dbl> 0.6666, 0.0000, 0.0000, 0.0000, 0.00…
## $ sca.spin <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.tric <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sca.spp <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sig.arge <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sig.cana <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sig.cora <dbl> 0.0000, 1.3332, 0.0000, 0.0000, 0.00…
## $ sig.doli <dbl> 5.3328, 3.9996, 1.3332, 2.6664, 6.66…
## $ sig.fusc <dbl> 0.0000, 0.0000, 11.9988, 6.6660, 0.0…
## $ sig.javu <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sig.line <dbl> 0.0000, 0.0000, 0.0000, 1.3332, 0.00…
## $ sig.puel <dbl> 0.0000, 0.0000, 0.6666, 0.0000, 0.00…
## $ sig.ptus <dbl> 0.0000, 0.0000, 0.0000, 1.9998, 0.00…
## $ sig..ptiss <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sig.stell <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ sig.vulp <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ mon.arge <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.auri <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ cha.afas <dbl> 63.3270, 89.9910, 39.3294, 57.3276, …
## $ cha.baro <dbl> 0.0000, 1.3332, 0.0000, 0.0000, 0.00…
## $ cha.benn <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.ephi <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ cha.citr <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.flav <dbl> 0.0000, 0.0000, 0.0000, 0.6666, 0.00…
## $ cha.line <dbl> 5.9994, 5.3328, 0.0000, 1.3332, 5.33…
## $ cha.lunu <dbl> 1.3332, 0.0000, 0.6666, 1.3332, 0.00…
## $ cha.ttus <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.mela <dbl> 4.6662, 1.3332, 1.9998, 0.0000, 0.66…
## $ cha.ocel <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.orna <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.pleb <dbl> 0.0000, 0.0000, 0.6666, 0.6666, 0.00…
## $ cha.rain <dbl> 3.3330, 13.9986, 1.3332, 3.3330, 3.9…
## $ cha.raff <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.seme <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cha.spec <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ cha.tlis <dbl> 1.3332, 2.6664, 0.0000, 0.0000, 0.00…
## $ cha.ulie <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ cha.vaga <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chm.rost <dbl> 7.3326, 3.3330, 1.9998, 1.3332, 6.66…
## $ cor.alti <dbl> 3.9996, 1.9998, 0.0000, 0.6666, 1.99…
## $ cor.chry <dbl> 0.0000, 1.3332, 0.0000, 0.6666, 1.99…
## $ hen.acum <dbl> 0.0000, 1.3332, 0.0000, 0.0000, 0.00…
## $ hen.mono <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ hen.vari <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ mct.stri <dbl> 0.0000, 0.0000, 0.0000, 0.0000, 0.00…
## $ prc.ocel <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cen.bico <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cen.bisp <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cen.nox <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cen.tibi <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ cen.vrol <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chd.doub <dbl> 5.9994, 3.9996, 0.6666, 1.3332, 7.99…
## $ chd.mere <dbl> 5.3328, 5.3328, 1.9998, 2.6664, 4.66…
## $ pmc.impe <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pmc.semi <dbl> 0.0000, 3.3330, 0.0000, 0.0000, 0.66…
## $ pmc.sext <dbl> 1.3332, 1.3332, 0.0000, 0.6666, 0.66…
## $ pmc.xant <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pyg.diac <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ abu.spp <int> 406, 48, 0, 0, 238, 0, 30, 112, 38, …
## $ acn.poly <int> 90, 48, 74, 58, 124, 22, 152, 78, 62…
## $ amb.aure <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ amb.cura <int> 0, 0, 4, 2, 0, 0, 0, 0, 0, 2, 0, 0, …
## $ amb.leuc <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ amp.spp <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ amp.peri <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pre.spp <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chl.labi <int> 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.ambo <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.alis <int> 60, 0, 12, 0, 0, 0, 138, 34, 16, 64,…
## $ chr.apes <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.marg <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.niti <int> 8920, 9680, 7340, 15840, 12210, 240,…
## $ chr.retr <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.tern <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.webe <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chr.xant <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chy.roll <int> 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chy.talb <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ chy.flav <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ das.spp <int> 12, 0, 0, 0, 0, 0, 0, 0, 0, 24, 4, 0…
## $ das.reti <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ das.trim <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ dis.spp <int> 0, 0, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ hgy.plag <int> 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0,…
## $ neg.mela <int> 14, 6, 2, 4, 12, 0, 4, 0, 8, 0, 0, 2…
## $ parma <int> 10, 14, 0, 0, 4, 2, 2, 2, 4, 0, 4, 0…
## $ neg.nigr <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pgy.dick <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pgy.lacr <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.adel <int> 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, …
## $ pom.ambo <int> 0, 0, 8, 0, 0, 0, 18, 6, 4, 32, 0, 0…
## $ pom.aust <int> 58, 140, 14, 54, 28, 0, 216, 156, 21…
## $ pom.bank <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.brac <int> 84, 0, 4, 6, 52, 0, 90, 14, 10, 0, 6…
## $ pom.chry <int> 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, …
## $ pom.coel <int> 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.gram <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.lepi <int> 8, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.molu <int> 602, 294, 426, 634, 542, 190, 564, 1…
## $ pom.naga <int> 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.phil <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.reid <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.trip <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.vaiu <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ pom.ward <int> 258, 224, 206, 338, 200, 456, 136, 1…
## $ ste.apic <int> 2, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, …
## $ ste.fasc <int> 8, 8, 4, 10, 8, 12, 0, 4, 8, 0, 4, 6…
## $ ste.livi <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ ste.nigr <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ oxy.long <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Anampses <int> 0, 6, 0, 2, 8, 0, 0, 0, 0, 0, 6, 2, …
## $ Coris <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Halichoeres <int> 14, 10, 20, 30, 8, 8, 26, 6, 10, 48,…
## $ Labroides <int> 18, 10, 10, 12, 16, 4, 10, 10, 4, 2,…
## $ Labropsis <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Psuedolabrus.guentheri <int> 20, 58, 20, 24, 26, 104, 16, 22, 20,…
## $ Stethojoulis <int> 12, 4, 16, 28, 2, 30, 16, 2, 0, 84, …
## $ Thalasoma <int> 24, 8, 12, 12, 24, 2, 30, 10, 16, 18…
## $ Macropharyngodon <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Labrichthys <int> 8, 0, 0, 6, 0, 2, 0, 2, 4, 4, 4, 0, …
## $ Total.fish.density <dbl> 11381.9272, 10805.9752, 8401.3122, 1…
## $ Total.fish.species.richness <dbl> 35.2, 29.4, 27.6, 26.4, 29.2, 17.8, …
## $ Prey.density <int> 10654, 10558, 8190, 17060, 13502, 10…
## $ Prey.biomass <dbl> 44.760041, 30.942239, 23.789262, 46.…
## $ Plectropomus.total.density <dbl> 5.9994, 3.9996, 17.3316, 7.3326, 10.…
## $ Plectropomus.total.biomass <dbl> 3.980542, 5.371889, 9.546584, 8.9263…
## $ Plectropomus.legal.density <dbl> 1.9998, 2.6664, 4.6662, 3.9996, 1.99…
## $ Plectropomus.legal.biomass <dbl> 2.4194697, 4.5827655, 5.0227904, 7.1…
## $ BE <dbl> 121.9966, 119.9976, 114.6630, 129.33…
## $ GRAZ <dbl> 72.6594, 47.9952, 36.6630, 26.6640, …
## $ COR <dbl> 103.3246, 119.3214, 45.9954, 71.9934…
## $ OM <dbl> 499.9996, 99.9996, 74.0000, 59.9998,…
## $ PL <dbl> 9000, 9680, 7366, 15842, 12210, 240,…
## $ CA <dbl> 521.9478, 45.9954, 65.9934, 43.9956,…
## $ PI <dbl> 7.9992, 6.6660, 25.9974, 10.6656, 12…
## $ FA <int> 290, 246, 218, 348, 212, 474, 138, 1…
## $ slope <dbl> 2.60, 2.64, 2.00, 2.40, 1.96, 2.20, …
## $ rugosity <dbl> 2.68, 3.08, 3.00, 3.00, 3.00, 2.68, …
## $ SCI <dbl> 7.040, 8.088, 6.000, 7.200, 5.880, 5…
## $ LCC_. <dbl> 60.0, 78.8, 77.6, 77.6, 57.2, 13.6, …
## $ LHC_. <dbl> 56.4, 68.8, 77.6, 75.6, 55.6, 13.6, …
## $ SC_. <dbl> 3.6, 10.0, 0.0, 2.0, 1.6, 0.0, 0.0, …
## $ MAC_. <dbl> 2.4, 0.4, 9.6, 4.4, 3.6, 4.0, 29.2, …
## $ SPO_. <dbl> 1.6, 0.0, 0.0, 0.4, 0.0, 0.0, 0.0, 0…
## $ Turf_. <dbl> 26.0, 17.6, 4.4, 11.6, 28.0, 71.2, 6…
## $ Unconsolidated_. <dbl> 9.6, 3.2, 8.4, 6.0, 11.2, 11.2, 7.2,…
## $ Benthic.richness <dbl> 6.6, 5.8, 2.2, 6.0, 3.8, 2.0, 3.4, 5…
## $ Coral_Morph.Diversity <dbl> 5.4, 5.6, 1.2, 5.0, 3.4, 1.4, 2.4, 4…
## $ ChlA <dbl> 0.45221510, 0.45221510, 0.41940367, …
## $ kd490 <dbl> 0.06519040, 0.06519040, 0.06261333, …
## $ SSTmean <fct> 23.80386083, 23.80357108, 23.7912496…
## $ SSTanom <dbl> -0.3130472, -0.3140394, -0.3159858, …
## $ wave.exposure.index <dbl> 0.007256172, 0.160425811, 0.00099007…
## $ Corrected.depth <dbl> 4.03, 6.39, 5.12, 5.47, 3.73, 3.02, …
## $ maxDHW <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, …
## $ deltaLHC <dbl> 14.0, 7.6, 6.4, 16.8, 8.0, 2.8, 4.8,…
## $ deltaMA <dbl> -6.8, -1.2, -3.6, 0.4, -6.4, -41.2, …
## $ Cyclone <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ LT.Fprimary <dbl> 0.087, 0.087, 0.087, 0.031, 0.025, 0…
## $ Exposure.to.primary.weeks <int> 2, 2, 2, 2, 3, 2, 3, 1, 1, 0, 3, 2, …
var.lookup = rbind(
data.frame(pretty.name='Total density', Field.name='Total.fish.density', Abbreviation='TFD', Family='gaussian', Type='Response', Transform='log', Groupby=''),
data.frame(pretty.name='Total species richness', Field.name='Total.fish.species.richness', Abbreviation='TFSR', Family='gaussian', Type='Response', Transform='I', Groupby=''),
data.frame(pretty.name='Benthic invertivores', Field.name='BE', Abbreviation='BE', Family='gaussian', Type='Response', Transform='log', Groupby=''),
data.frame(pretty.name='Grazers', Field.name='GRAZ', Abbreviation='GRAZ', Family='gaussian', Type='Response', Transform='log', Groupby=''),
data.frame(pretty.name='Corallivores', Field.name='COR', Abbreviation='COR', Family='gaussian', Type='Response', Transform='log', Groupby=''),
data.frame(pretty.name='Omnivores', Field.name='OM', Abbreviation='OM', Family='gaussian', Type='Response', Transform='log', Groupby=''),
data.frame(pretty.name='Planktivores', Field.name='PL', Abbreviation='PL', Family='gaussian', Type='Response', Transform='log', Groupby=''),
data.frame(pretty.name='Carnivores', Field.name='CA', Abbreviation='CA', Family='gaussian', Type='Response', Transform='log', Groupby=''),
data.frame(pretty.name='Piscivores', Field.name='PI', Abbreviation='PI', Family='gaussian', Type='Response', Transform='log', Groupby=''),
data.frame(pretty.name='Farmers', Field.name='FA', Abbreviation='FA', Family='gaussian', Type='Response', Transform='log', Groupby=''),
data.frame(pretty.name='Plectropomus total density', Field.name='Plectropomus.total.density', Abbreviation='PTD', Family='gaussian', Type='Response', Transform='log', Groupby=''),
data.frame(pretty.name='Plectropomus total biomass', Field.name='Plectropomus.total.biomass', Abbreviation='PTB', Family='gaussian', Type='Response', Transform='log', Groupby=''),
data.frame(pretty.name='Plectropomus legal density', Field.name='Plectropomus.legal.density', Abbreviation='PLD', Family='gaussian', Type='Response', Transform='log', Groupby=''),
data.frame(pretty.name='Plectropomus legal biomass', Field.name='Plectropomus.legal.biomass', Abbreviation='PLB', Family='gaussian', Type='Response', Transform='log', Groupby=''),
data.frame(pretty.name='Region', Field.name='REGION', Abbreviation='REGION', Family=NA, Type='Predictor', Transform='I', Groupby=''),
data.frame(pretty.name='NTR Pooled', Field.name='NTR.Pooled', Abbreviation='NTR.Pooled', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='LHC % (live hard coral cover)', Field.name='LHC_.', Abbreviation='LHC', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='SC % (soft coral)', Field.name='SC_.', Abbreviation='SC', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='MA % (macroalgal cover)', Field.name='MAC_.', Abbreviation='MA', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Turf %', Field.name='Turf_.', Abbreviation='TURF', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Unconsolidated %', Field.name='Unconsolidated_.', Abbreviation='UC', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Benthic richness', Field.name='Benthic.richness', Abbreviation='BR', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Coral morphological diversity', Field.name='Coral_Morph.Diversity', Abbreviation='CMD', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Slope', Field.name='slope', Abbreviation='SLOPE', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Rugosity', Field.name='rugosity', Abbreviation='RUG', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='SCI (structural complexity)', Field.name='SCI', Abbreviation='SCI', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Chla', Field.name='ChlA', Abbreviation='CHL', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Kd490', Field.name='kd490', Abbreviation='KD490', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='SST mean', Field.name='SSTmean', Abbreviation='SSTMEAN', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='SST anom', Field.name='SSTanom', Abbreviation='SSTANOM', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Wave exposure', Field.name='wave.exposure.index', Abbreviation='WAVE', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Corrected depth', Field.name='Corrected.depth', Abbreviation='DEPTH', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Max DHW', Field.name='maxDHW', Abbreviation='DHW', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Cyclone', Field.name='Cyclone', Abbreviation='CYCLONE', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Exposure to primary weeks', Field.name='Exposure.to.primary.weeks', Abbreviation='EXP', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Prey density', Field.name='Prey.density', Abbreviation='PREY.DENSITY', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Prey biomass', Field.name='Prey.biomass', Abbreviation='PREY.BIOMASS', Family=NA, Type='Predictor', Transform='I', Groupby='Region')
)
save(var.lookup, file='data/var.lookup.RData')
names = with(var.lookup, setNames(as.character(Field.name), Abbreviation))
## exclude those whose names equal their values otherwise there will be duplicate fields created in the fish data
names = names[names(names)!=names]
## Create duplicates of the fields that are not already abbreviated
fish = fish %>%
mutate(SSTmean=ifelse(as.character(SSTmean)=='#N/A',NA,as.numeric(as.character(SSTmean)))) %>%
mutate_at(as.character(var.lookup$Field.name), list(A=~I)) %>%
rename(!!! gsub('(.*)','\\1_A',names)) %>%
mutate(REGION=factor(REGION, levels=c('Palm','Magnetic','Whitsunday','Keppel')))
save(fish, file='data/fish.RData')
Species lists
For each model, Dani has compiled a list of appropriate Species
speciesList = read.csv('data/Species selection.csv', strip.white=TRUE)
speciesList %>% glimpse
## Observations: 248
## Variables: 6
## $ Palms.model <fct> sca.rivu, sig.doli, sca.flav, chs.sord, sca.gh…
## $ Magnetic.model <fct> sca.rivu, sca.ghob, aca.duss, aca.bloc, sig.do…
## $ Whitsunday.model <fct> sca.rivu, sig.doli, sca.flav, sca.nigr, chs.so…
## $ Keppel.model <fct> sca.rivu, sca.ghob, sig.arge, sig.fusc, sig.do…
## $ X <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ All.regions.model <fct> aca.line, aca.ncus, nas.brach, nas.litu, cal.c…
speciesList = speciesList %>% dplyr::select(-X) %>% dplyr::rename(`all.model`=`All.regions.model`, `Palm.model`=`Palms.model`)
For consistency with the functional groups analyses, I will turn this list into a list of Species that indicates which location models the species should be involved in.
speciesNames = levels(unlist(speciesList))
speciesNames = speciesNames[speciesNames!=""]
## The names dont quite line up between the two data sets
## The following lines are to help them match
speciesNames = gsub(' ','.',speciesNames)
speciesNames = speciesNames[speciesNames!="pse.tuka"]
speciesList = sapply(speciesNames, function(x) colSums(speciesList==x[1]), simplify=FALSE,USE.NAMES = TRUE)
save(speciesList, file='data/speciesList.RData')
var.lookup.species = do.call('rbind',
lapply(speciesNames, function(x) data.frame(pretty.name=x, Field.name=x,
Abbreviation=x, Family='poisson',
Type='Response', Transform='I',
Groupby='')
)
)
var.lookup.species = rbind(
var.lookup.species,
data.frame(pretty.name='Region', Field.name='REGION', Abbreviation='REGION', Family=NA, Type='Predictor', Transform='I', Groupby=''),
data.frame(pretty.name='NTR Pooled', Field.name='NTR.Pooled', Abbreviation='NTR.Pooled', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='LHC % (live hard coral cover)', Field.name='LHC_.', Abbreviation='LHC', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='SC % (soft coral)', Field.name='SC_.', Abbreviation='SC', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='MA % (macroalgal cover)', Field.name='MAC_.', Abbreviation='MA', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Turf %', Field.name='Turf_.', Abbreviation='TURF', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Unconsolidated %', Field.name='Unconsolidated_.', Abbreviation='UC', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Benthic richness', Field.name='Benthic.richness', Abbreviation='BR', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Coral morphological diversity', Field.name='Coral_Morph.Diversity', Abbreviation='CMD', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Slope', Field.name='slope', Abbreviation='SLOPE', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Rugosity', Field.name='rugosity', Abbreviation='RUG', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='SCI (structural complexity)', Field.name='SCI', Abbreviation='SCI', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Chla', Field.name='ChlA', Abbreviation='CHL', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Kd490', Field.name='kd490', Abbreviation='KD490', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='SST mean', Field.name='SSTmean', Abbreviation='SSTMEAN', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='SST anom', Field.name='SSTanom', Abbreviation='SSTANOM', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Wave exposure', Field.name='wave.exposure.index', Abbreviation='WAVE', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Corrected depth', Field.name='Corrected.depth', Abbreviation='DEPTH', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Max DHW', Field.name='maxDHW', Abbreviation='DHW', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Cyclone', Field.name='Cyclone', Abbreviation='CYCLONE', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Exposure to primary weeks', Field.name='Exposure.to.primary.weeks', Abbreviation='EXP', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Prey density', Field.name='Prey.density', Abbreviation='PREY.DENSITY', Family=NA, Type='Predictor', Transform='I', Groupby='Region'),
data.frame(pretty.name='Prey biomass', Field.name='Prey.biomass', Abbreviation='PREY.BIOMASS', Family=NA, Type='Predictor', Transform='I', Groupby='Region')
)
save(var.lookup.species, file='data/var.lookup.species.RData')
names.species = with(var.lookup.species, setNames(as.character(Field.name), Abbreviation))
## exclude those whose names equal their values otherwise there will be duplicate fields created in the fish data
names.species = names.species[names(names.species)!=names.species]
## Create duplicates of the fields that are not already abbreviated
fish = read.csv('data/Fish benthic physical sitelevel 2019.csv', strip.white=TRUE)
fish.species = fish %>%
mutate(SSTmean=ifelse(as.character(SSTmean)=='#N/A',NA,as.numeric(as.character(SSTmean)))) %>%
mutate_at(as.character(var.lookup.species$Field.name), list(A=~I)) %>%
rename(!!! gsub('(.*)','\\1_A',names.species)) %>%
mutate(REGION=factor(REGION, levels=c('Palm','Magnetic','Whitsunday','Keppel')))
save(fish.species, file='data/fish.species.RData')
load('data/var.lookup.species.RData')
resp.lookup = var.lookup.species %>% filter(Type=='Response') %>% droplevels
for (s in speciesNames) {
#print(s)
resp=s
cat(paste('## ',resp.lookup$pretty.name[resp.lookup$Abbreviation==resp],' {.tabset .tabset-pills} \n\n'))
EDA_histograms(var=s, dat=fish, var.lookup=var.lookup.species)
cat('\n\n')
EDA_density(var=s, group='REGION', dat=fish, var.lookup=var.lookup.species)
cat('\n\n')
}
It is clear that all species counts have a long right tail and a large number of zeros. There are no zero-inflated regression tree models as such. If necessary, we could spilt the analyses up into Presence/Absence (binary) and Poisson (>0) models.
For exploratory data analyses of explanatory (predictor) variables, see here.
formulas = list(
all = ~REGION + NTR.Pooled + LHC + SC + MA + TURF + UC + BR + CMD + SLOPE + RUG + SCI + CHL + KD490 + SSTMEAN + SSTANOM + WAVE + DEPTH + DHW + CYCLONE + EXP,
Palm = ~ NTR.Pooled + LHC + SC + MA + TURF + UC + BR + CMD + SLOPE + RUG + SCI + CHL + KD490 + SSTMEAN + SSTANOM + WAVE + DEPTH + DHW + CYCLONE + EXP,
Magnetic = ~ NTR.Pooled + LHC + SC + MA + TURF + UC + BR + CMD + SLOPE + RUG + SCI + CHL + KD490 + SSTMEAN + SSTANOM + WAVE + DEPTH + DHW + CYCLONE + EXP,
Whitsunday = ~ NTR.Pooled + LHC + SC + MA + TURF + UC + BR + CMD + SLOPE + RUG + SCI + CHL + KD490 + SSTMEAN + SSTANOM + WAVE + DEPTH + DHW + CYCLONE + EXP,
Keppel = ~ NTR.Pooled + LHC + SC + MA + TURF + UC + BR + CMD + SLOPE + RUG + SCI + CHL + KD490 + SSTMEAN + SSTANOM + WAVE + DEPTH + DHW + CYCLONE + EXP
)
load('data/fish.species.RData')
load('data/var.lookup.species.RData')
resp.lookup.species = var.lookup.species %>% filter(Type=='Response') %>% droplevels
pred.lookup.species = var.lookup.species %>% filter(Type=='Predictor') %>% droplevels
groupings.all.species = vector('list', nrow(pred.lookup.species))
names(groupings.all.species) = pred.lookup.species$Abbreviation
groupings.region.species = vector('list', nrow(pred.lookup.species))
names(groupings.region.species) = pred.lookup.species$Abbreviation
for (i in 1:nrow(pred.lookup.species)) {
pred = as.character(pred.lookup.species[i,'Abbreviation'])
groupings.all.species[[pred]] = ifelse(pred=='REGION',NA,'REGION')
groupings.region.species[[pred]] = ifelse(pred=='NTR.Pooled',NA,'NTR.Pooled')
}
groupings.all.species = do.call('c',groupings.all.species)
groupings.region.species = do.call('c',groupings.region.species)
gfun = function(f, form) {
#print(f)
#print(form[[f]])
form = attr(terms(form[[f]]),'term.labels')
if (f=='all') return(groupings.all.species[form])
else return(groupings.region.species[form])
}
analyses.species = vector('list',nrow(resp.lookup.species))
names(analyses.species) = resp.lookup.species$Abbreviation
for (i in 1:nrow(resp.lookup.species)) {
resp = as.character(resp.lookup.species[i,'Abbreviation'])
fun = as.character(resp.lookup.species[i,'Transform'])
fam = as.character(resp.lookup.species[i,'Family'])
analyses.species[[resp]] = list('formulas' = lapply(formulas, function(f) update(f, paste0(fun,'(',resp,') ~.'))),
'family' = fam)
#if (resp %in% c('PI','PTD','PLD')) analyses[[resp]][['formulas']] = lapply(analyses[[resp]][['formulas']], function(f) f=update(f, .~.+PREY.DENSITY))
#if (resp %in% c('PTB','PLB')) analyses[[resp]][['formulas']] = lapply(analyses[[resp]][['formulas']], function(f) f=update(f, .~.+PREY.BIOMASS))
analyses.species[[resp]][['groups']] = sapply(c('all','Palm','Magnetic','Whitsunday','Keppel'), gfun, form=analyses.species[[resp]][['formulas']], USE.NAMES = TRUE,simplify=FALSE)
}
save(analyses.species, file='data/analyses.species.RData')
library(ssh)
session = ssh_connect('mlogan@hpc-login', keyfile='~/.ssh/hpc_id_rsa')
scp_upload(session,
files=c('FISH_functions.R',
'FISH_HPC_gbm.R',
'FISH_HPC_gbm.batch'),
to = "~/tmp/fish", verbose = TRUE)
scp_upload(session,
files=c('data/analyses.species.RData',
'data/fish.species.RData',
'data/speciesList.RData',
'data/var.lookup.species.RData'),
to = "~/tmp/fish/data", verbose = TRUE)
ssh_exec_wait(session,
command = "qsub -l mem=1gb -l nodes=1:ppn=20 ~/tmp/fish/FISH_HPC_gbm.batch\n")
ssh_exec_wait(session, command = "qstat -a")
#ssh_exec_wait(session, command = "qdel 71728")
ssh_disconnect(session)
session = ssh_connect('mlogan@hpc-login', keyfile='~/.ssh/hpc_id_rsa')
scp_download(session,
files=c('/export/home/l-p/mlogan/tmp/fish/output/figures/*_in.p*'),
to = "output/figures")
scp_download(session,
files=c('/export/home/l-p/mlogan/tmp/fish/data/pred.1*.RData'),
to = "data/")
scp_download(session,
files=c('/export/home/l-p/mlogan/tmp/fish/data/stats*.RData'),
to = "data/")
scp_download(session,
files=c('/export/home/l-p/mlogan/tmp/fish/data/thresholds*.RData'),
to = "data/")
# Note, for some reason the above returns an error, even through it is successful
ssh_disconnect(session)